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1 . Introduction 


In the present Report, I discuss a technique for the analysis 
of two-dimensional, inviscid, compressible, transonic flows. Em- 
phasis is put on the treatment of shocks. I will support my argu- 
ments with some numerical experiments on two-dimensional flows in 
ducts . 

Whatever the range of Mach numbers, indeed, such flows pro- 
vide a fertile ground for testing of numerical techniques and for 
the analysis of not-so-well-known features of steady and unsteady 
gasdynamical patterns. Topics of theoretical interest are: For- 

mation and evolution of shocks, regular and Mach reflections, sta- 
bility of shock patterns and (not considered in the present Re- 
port) interaction of unsteady internal shock patterns with the 
impinging flow. From a numerical standpoint, the relative merits 
and disadvantages of flux-difference splitting methods and of the 
A-scheme with shock tracking can be tested. 

The numerical scheme must: 

a) be robust in unsteady, transonic regions, 

b) properly represent the chosen inlet and outlet boundary 
models , 

c) properly satisfy the boundary conditions on rigid walls, 

d) provide accurate predictions of shock formation, regard- 
less of the number and location of shocks, 

e) provide an accurate analysis of shocks, their motion and 
interactions. 

It is desired to satisfy the above requirements with a code, 
simple to understand, easily extendable to other problems, and as 
fast as possible. For a time-dependent code, where the advance in 
time is controlled by the CFL rule, the latter condition requires 
the use of coarse meshes, the reduction of calculations at any 
point to what is strictly necessary (elimination of redundancies), 
the calculation itself streamlined to a minimal number of opera- 
tions (particularly multiplications and divisions), and the ab- 
sence of logical or arithmetic branches ( IF-statements) . 


In the present paper, I intend to outline a scheme which 



seems to satisfy the above conditions. 


2. The problem 

Let an inlet be schematized as a duct, running from x=-x 0 to 
x=x c . The shape of the lower wall is defined by y=b(x). The 
shape of the upper wall is defined by y=c(x). To simplify the 
analysis, we exclude any possibility of transonic interaction of 
the inlet flow with an infinite external flow; consequently, we 
can use a simple model to represent the external flow at the in- 
flow boundary of the duct. In the same spirit of simplicity, we 
assume that the inlet itself is preceded by a straight channel and 
followed by another straight channel (Fig. 1). 



In the numerical analysis, the flow will be computed between 
two sections, located at x=-Xj and x=x 2 . The flow problem can be 
modeled in three different ways: 

1) The duct is permanently at rest. The gas in the duct is 
initially at rest. A diaphragm is abruptly broken at x=x 2 , at 
time t=0, letting the gas in the duct communicate with an infinite 
cavity where the pressure is prescribed and constant. If the 
latter is lower then the initial pressure of the gas in the duct, 
a flow is generated and it will eventually reach a steady state. 

2) The duct is permanently at rest. The gas in the duct is 
initially at rest. A normal shock wave reaches the duct from the 
left. We can think of it as generated by the rupture of a 
diaphragm located far to the left of -Xj, or by the push of a pis- 
ton, started from rest with a constant speed at a time t<0. In 
either case, the shock is followed by a constant state. The shock 
itself travels from left to right within the duct and eventually 
disappears to the right of x 2 . 

3) The duct is initially at rest. The gas in the duct is 
also at rest. At t=0, the duct is impulsively set into motion 
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from right to left, in a direction parallel to the x-axis, with a 
velocity which is then maintained constant. The gas flow will be 
described by an observer moving with the duct, who will use the 
same frame of reference shown in Fig. 1. This, somewhat more ar- 
tificial, way of starting the calculation may be interpreted as 
the limit of the more physical procedure of accelerating the duct 
from rest, the entire accelerating process taking place in one 
single computational step. 


3. Equations of motion 

As usual, we assume a pressure, a density and a length as 
reference values (P ref » p f » x ref * respectively) . The reference 

velocity, temperature and time are thus: 

u = Hp T/p 7 , 0 = p ,/(Rp , ), t = x ,/u . (1) 

ref >\| ref ref ref ref ref ref ref ref 

where R is the gas constant. A nondimensional entropy, s, is de- 
fined as 


s = Y(y=T) [ln p ~ Y ln 

Consequently, 


( 2 ) 


0 = exp[ (Y-1 ) (P/Y + s)] 
a 2 = Y0 

where Y is the ratio of specific heats. 

The equations of motion can be written in the form: 
P fc + q. VP + YV.q = 0 
q t + (q.V)q + —■ VP = ° 
s t + q.Vs = 0 


(3) 

(4) 


(5) 


where P = ln p, q = Ui + Vj , and i, j are the unit vectors of the 
x and y axes, respectively. From (3) and (4) it follows that 


Y da v , 
dP = — — - Y ds 
6 a 


( 6 ) 


where 


6 = 


Y-1 


(7) 
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Therefore, (5) may be replaced by 


t a 

— + q.V -r + a V.q - a s - a q.Vs = 0 
0 c z 

q. + (q.V)q + a V 7 - a 2 Vs = 0 ( 8 ) 

z 0 

s + q.Vs = 0 

which are the obvious extension of Eqs. (4) of Ref. 1 to non- 
isentropic flow. 

One can use the same procedure of Ref. 1, defining two arbi- 
trary, but orthogonal, unit vectors n and \ at any point, and cal- 
ling a 0 the angle formed by n and i. Using the same notations as 
in Ref. 1, in particular 

g = q.Va 0 (9) 

F = a k>q. Va c ( 1 0 ) 

where k = i>j, and also 

C = cos a 0 , S = sin a e (11) 

a 

P!=y+n.q, Ai=q+an 

Ps = T " n.q, A £ = q - an 

a (12) 

p 3 = j + i *Q» A 3 = q + a-i 

a 

Pl * = i " T,q ’ Ai, = q - ai 
the equations of motion take on the form: 

■7 a + l A. . ( Vp . - a Vs) - 2q.V ~ + 4aq.Vs + 2F = 0 

0 v 1 i 0 

2n.a + Aj.(Vpj- a Vs) - A 2 .(Vp a - a Vs) - 2£q.i = 0 

X (13) 

2i.q t + A 3 . (Vp 3 - a Vs) - A^.fVp,,- a Vs) + 2gq.n = 0 

s fc + q.Vs = 0 

We will also introduce the velocity components: 

u = CU + SV 

(14) 

v = - SU + CV 

along the n- and t-directions, respectively. 
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3.1. Computational grid 


We will choose a computational grid by defining two computa- 
tional variables, X and Y, as follows: 


X = x 

y-b(x) 
e(x )-b( x) 


( 15 ) 


from which: 

Y = — — , Y = Y [b'(Y-1)-c'Y] (16) 

y c-b x y 

(primes denote differentiations with respect to x). The grid will 
result as shown in Fig. 2. It is obviously not orthogonal between 
-x 0 and x 0 . We will assume n to be the direction of the Y = con- 
stant lines in the grid. Consequently, 

C = Y /v , S = - Y/v , v = (Y 2 +Y 2 ) (17) 

y ^ y 

a = - C 2 [b”(Y-1 )-c"Y-Y (c'-b')] 
ox x (18) 

a = v 2 /Y (c'-b') 

oy y 

(5 = Ua + Va (19) 

ox oy 

F = a (Ua - Va ) (20) 

oy ox 

Also, note that 

CY - SY = v , CY + SY = 0 (21) 

y x x y 


Therefore, 


VY + UY = w 

y x 


( 22 ) 


-Xa 
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3.2. Generalized Riemann variables 


In the spirit of Ref. 1, we introduce generalized Riemann 
variables in the form: 

R? - | - u. kJ - f - u, R* - f . v, R, - | - V 
and corresponding expressions for eight new coefficients: 

A* = U+Ca, A* = U-Ca, A 3 = U, aJ = U-Sa, A? = U+Sa 
A^ = ( v+a)Y /C, X Y Z = (v-a)Y /C, A 3 = vY /C 

J J ■ J 

Finally, the following terms are defined: 


f* - - »X x -as x >. 

f X ,X,_Y „ » 

f\ - ~ X^(R^x ^ s x J 

fj = “ _as y^ ’ 


,x,_x s ,x ,x a x . 

A 2 R 2X aS X ’ fs " 3 { as x } 


X 

A 5 (R 2X 3 S X^ ’ 

.X ,x 

fY = - 

f 6 - A 3 s x 

►“h 

ro F-< 
1! 

Y Y 

A 2 (R 2Y _aS Y^’ 

,Y , Y 

f 3 = " * 3 ^ 

-Y 

, Y 


f «♦ ” ” 

A 3 S Y 



fs = 2 gv, f 3 = - 2 gu 


A new form of the equations of motion results: 

,X _Y 
s t = f 6 +f 4 

6 , ,X ,X ,Y ,Y X X X Lx 
a t = 2 (f 1 +f 2 +f 1 +f 2 +f 4 +f 5 2 f 3 +f 1 } 
1 ' X ^X 0 ^Y 
U t 2 f 1 ~ f 2 + 2 f 3 +f 2 

1 , Y Y .X XL. 

V t = 2 (f r f 2 +f 4 _f 5 f 3 


- 6 as, 


(23) 

(24) 

(25) 

(26) 


As pointed out in Ref. 1, in the new formulation there are 

x 

terms (denoted by f.) which contain only X-derivatives, terms 
Y 1 

(denoted by f. ) which contain only Y-derivatives , and local terms 

1 L 

(denoted by f\.). Such a splitting provides a clear definition of 

the points to be taken into account in the discretization of 
derivatives. According to the arguments exposed in Ref. 1 , each 
X- or Y-derivative is to be approximated using only two points. 
One is the point to be updated, the other is the point on the X- 
or Y-axis, respectively, near the point to be updated, on the side 

X X 

from which the line X=A^t or Y=A^t proceeds. Specifically, for 

the X-terms, the point to be used lies to the left of the point to 

X X 

be updated if A^>0, and to the right of it if AXO. For the Y- 

terms, the point to be used lies below the point to be updated if 

Y Y 

A.>0, and above it if A.<0. 

1 1 
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4. Two- level integration scheme 

To avoid unnecessary branching options (which are time 
consuming on the computer), we will extend the computational re- 
gion by one additional row of points all around the boundary. 
This will allow certain derivatives to be approximated by differ- 
ences taken from outside, without violating any physical principle 
since such terms are computed but never used, as explained in the 
next Section. 

At all points, the values of a, u, v and s are updated using 
a two-level scheme, formulated as follows: 


4.1. Predictor level . 

X y l 

All values of f^, f^, and f^ are computed, approximating any 
X-derivative according to the formula: 


where 


^X = 


J = n 


<j> . . 
J+1 »m 

^j .m 

AX 


n if 

A X <0 

l 

n-1 if 

A X >0 


and any Y-derivative according to the formula: 


where 


*n,j+1 



AY 


j = m 

if 

Y 

A . <0 

l 

j = m-1 

if 

A Y >0 


All coefficients are approximated accordingly as 


A 


X 


A 


Y 


1 r 

* tt L A . + A 

2 j ,m 

1 r J 

= -x [ A . + A 

2 n,j 


X 

j + 1 ,m 
Y 

n,j+1 


] 

] 


Note that all terms are defined according to their physical 
domains of dependence at all points. 


made 


X XX 

At nodes bracketing a transition from A 2 <0 to A 2 >0, A 2 is 

X 

equal to the local value, A . This device (which can be 

n ,m 
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applied without using IF-statements) prevents the formation of ex- 
pansion shocks. In the presence of shock, similar devices must be 
used in their vicinity. They will be explained in Section 6 . 3 . 

X Y L 

All values of f\ , f\. and f^ so def ined at the beginning of 

(k) 

the computational step will be denoted by f for brevity. They 
will be used to advance half a step using (23) ,(24) ,(25) and (26) 
in conjunction with 


k+1/2 

a 


k 

a + 


If (k) At/2 


(27) 


and similar formulas. 


4.2. Corrector level . 

All f's are recomputed using the values updated by (27) at 
half step. Then, in order to maintain the proper domain of depen- 
dence at the corrector level, we denote the new f’s by f’s and 
correct them using the formulas: 


X(k+1/2) ,X _ 1_ f x < k ) 

2 2j-n+1,m 

Y(k+1/2) ~Y J_ Y(k) 

2 n ,2j-m+1 

f L(k+1/2) _ f L _ _1_ L(k) 

! 2 n,m 


Finally, formulas similar to (27) are applied: 


a 


k+1 


k+1/2 v.ik+l/Z) 
a + /f 


At 


(28) 


X X 

A remark is in order as regards f„ and f 5 (typical contri- 
butions due to the non-orthogonality of the grid). They appear in 

two combinations only: g=f ,,^+f 5 X -2f 3 ^ and h=f,,^-f 5 ^. With U>0 and 
U>Sa (a most likely occurrence), all derivatives in g and h are 
approximated by upstream differences. at the corrector level, 
only g and h (both computed at n-1 ) are needed in the updating of 

X 

a and v, respectively. If U<Sa, f,, is approximated by downstream 
differences. We will make use of these remarks in Sections 5 and 

6 . 1 . 


5. Boundary conditions 

Boundary conditions must be applied at four boundaries: 1) 

lower wall, 2) upper wall, 3) entrance, and 4) exit. On the 

x 

walls, no special provision has to be made for the f terms or for 
L Y 

the f terms. As regards the f terms, note that v=0 on both 

Y 

lower and upper wall. Therefore, f 3 =0. All local terms also 


8 



Y 

vanish. On the upper wall, f, is computed using information from 

Y 

inside the flow; therefore, it is correct; f 2 , instead, must be 
defined using the condition, v^=0. On the lower wall, the oppo- 

Y 

site is true. Consequently, f 2 will be determined on the upper 
wall by the equation 

fa = f] * f\ - fs (29) 

Y 

which descends from (26), and fj will be determined on the lower 
wall by 

f^=f 2 (30) 

which also descends from (26), with the additional, obvious con- 

X X 

sideration that f,, -f 5 =0 (since the wall lies on the x-axis and 
S-0, A^=Aj, R*=R 2 ). 

At the entry and exit boundaries, the local terms vanish, be- 
cause the walls are parallel. On the entry boundary the flow may 
be supersonic or subsonic. If it is supersonic, all values can be 
prescribed, and all f-s vanish identically. If the flow is sub- 

Y 

sonic, we assume that v=0 and, consequently, f 3 =0 at the left 

boundary. The entropy is constant in the impinging flow (to the 

X Y XXX 

left of -X!); therefore, f 6 =f„=0. Moreover, A 3 =A 1( =A e ; consequent- 
ly, 

g = f l + f* “ 2f* =0, h = f* -f* =0 (31) 

The impinging flow has a clearly defined total temperature. It is 
convenient, thus, to impose that the total temperature is also 
constant at x=-x } all across the boundary. From the definition of 
total speed of sound, 

2 2 2 

af = a + 6q (32) 

written in the form: 

a + 6 q q fc = 0 (33) 

and the assumption that v=v t =0 at x=-Xj , (24) and (25) with the 
first of (31) yield: 

(a+u)f* + (a-u)!^ + a(f^+ f*) = 0 (34) 

X 

from which f i can be evaluated. Note that any value can be used 
X XX 

for f^, provided that f 3 and f 3 are defined to satisfy (31). Such 
a procedure is correct because only the values of g and h computed 
at the boundary are used, at the corrector level, to correct the 
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predicted values on the line next to the boundary. 

On the right boundary, if the flow is supersonic, no special 
provision is to be taken (all X-derivatives are automatically ap- 
proximated by upstream differences). If the flow is subsonic, we 
will simply impose that the pressure, p , is constant in time. 
Using (6), the condition is found: 

a /6 - a = 0 (35) 

Note that g=0 at the exit, where S=0. Therefore, we can apply 
(35) after setting g equal to zero, and, from (35), we find the 
condi tion, 

f X + f 2 + f*+ t\ = ° (36) 

x 

which allows f £ to be evaluated from known quantities. 


5.1. Resetting . 

y 

Terms, such as f 2 in (30), re-evaluated at the boundaries, 
are used in (27) and (28) to update the boundary points them- 
selves. In principle, v, originally equal to zero on the walls, 
should remain equal to zero because (30) assures the vanishing of 
v. . Similarly, a c and p gx should remain constant in force of (3*0 

and (35). In practice, it may not be so because the updating of 
v, a c and p gx is affected by almost imperceptible truncation er- 
rors in time. After a number of steps, however, one can observe a 
departure from the original values, producing an increase or de- 
crease in total energy and/or a non-vanishing v (expressing an ad- 
dition or loss of mass through the wall). It is necessary, there- 
fore, to reset certain quantities to maintain a c , p gx and v con- 
stant at the entry, exit and wall boundary, respectively. The 
task is accomplished, at x=-Xj , by computing 

R 2 X = a/6 - u (37) 

from the available values of a and u and using it, together with 
( 32 ) , to obtain: 

1 r rn X ,6+1 2 _, n X. 2.1/2-, ,_ D v 

u = [-6R* + (-£- a c 2 - 6 (R 2 ) 2 ) ] (38) 

Then, a follows from (37). At x=x 2 , 

a = /y exp [6(ln p/Y + s)] (39) 

x 

and u is recomputed from the definition of Rj : 

u = Rj X - a/6 (40) 

Finally, at the upper wall, 
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( 41 ) 


_ -A. 

Rj = a/c + v 


is recomputed from the available values of a and v and the 
tion, v=0 is used to obtain a new value for a: 


condi - 


a = 5 R j 


(42) 


U and V are obtained from (14), with v=0 and u accepted as comput- 
ed. 


6. Shocks 

The analysis of a flow containing shocks cannot be made using 
the technique described in Sections 3 through 5, as is. Indeed, 
the equations of motion written in the form (8) assume that the 
entropy is carried by the moving flow elements. Therefore, such 
equations cannot account for the rise in entropy occurring at 
shocks. Oblique shocks with low normal Mach numbers can be prop- 
erly captured by the technique . If normal shocks or shocks fol- 
lowed by a subsonic region occur, additional care must be taken to 
handle them correctly. Other techniques [2,3,4], seem to reach 
the goal of capturing all kinds of shocks. This is done by using 
at every point the same computational scheme which would be needed 
at shock points only. In the present approach, the technique for 
all points remains the same as described above. All shock points, 
however, must be tracked and the Rankine-Hugoniot conditions used 
to update the grid points downstream of a shock. 

The updating of a shock point is very simple and straightfor- 
ward, once the direction of the normal to the shock is known be- 
cause, along the normal, the one-dimensional technique [5] can be 
used. Finding the normal to a shock at any shock point, however, 
requires a tracking of all the shocks, which involves a certain 
amount of logic and bookkeeping. We will begin with discussing 
the updating at any shock point, assuming that the direction of 
the normal thereon is already known. 


€i 



Figure 3.- Shock wave passing through Y=constant line. 
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6.1. Updating a_ shock point . 

Let a (Fig. 3) be a shock wave, separating a low-pressure re- 
gion from a high-pressure region. Let A and B be two points, on a 
Y=constant line, on either side of the shock. We know their (com- 
mon) position, as well as the direction of the normal, defined by 
the unit vector, N. Both at the predictor and at the corrector 
level, the calculation of the shock follows the updating of all 
grid points, made following the procedure detailed in Section *T 
(with a few provisoes to be explained in Section 6.3). 

From now on, we follow. the analysis exposed in Ref. 5, Sec- 
tion 9 for normal shocks. Let Nj and N £ be the cosine and sine of 
the angle between N and n, and 

u = q . N = uNj + vNj (43) 

the component of the velocity normal to the shock; let us call v 
the tangential component. 


v = — uN £ + vNj (44) 

The parameter, Z, defined by (30) in Ref. 5, is redefined here as 
l = 6 (Rg B + u A ) /a A = + U/a A )[u] (45) 

where [] denotes a jump across the shock. Note that 


[u] = [u] Nj (46) 

since [u]=[uN a -vN £ ] from (43) and (44), but [v]=0. The normal re- 
lative Mach number of the shock is 


M = 



(47) 


where WN is the velocity of the shock as a locus of discontinui- 
ties. From the Rankine-Hugoniot conditions we obtain: 


1-M 

U B U A + a A ( 1 +6 )M 


a B S A 
1 


A 


(YM 2 -6) (1+6M 2 ) 


(1+6) |M 


ri YM z -6 v 1+6M 2 , 
S " S A + 2Y6 ^ ln 1+6 + T ln (l+ejM 2 - 1 


(48) 

(49) 

(50) 



(51) 


From the last term in (45), using (48) and (49), we obtain 
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V 


(52) 


L__ [ v|(7M*-0 (1+6M Z ) + 6(M Z -1)N,] 

( 1 +6 )M L ^ ' ' ’ ' 0 

From a given I, H is thus obtainable. Since (52) cannot be in- 
verted, it is convenient to proceed as follows. With 

g» = \|(16 y- 6)(1+166)]/[4(1+6)] , g £ = (156)/[4( 1+6)] (53) 

l c = Si + (54) 


a first guess for M is written in the form: 


M = Gs+GjL, Gj = 3/(Z 0 -1),G s = 1-G X (55) 

and the approximation is improved iteratively, until E gets sta- 
bilized to within a given tolerance (0.00001 on an IBM 370, 
0.000001 on a CRAY 1). Three to five iterations are sufficient. 
The shock velocity is then obtained from (47), and the Rankine- 
Hugoniot conditions can be applied to update the values at B. The 
abscissa of the shock point is obtained by adding 


to its current value. 


Ax, 


W At 


cosa 0 

N, 


(56) 


With different degrees of sophistication, we can assume that 
A and B coincide with the grid nodes bracketing the shock point 
(an acceptable approximation if the grid is very fine) or extrapo- 
late upstream values on the same Y=constant line to define values 
at A and downstream values to define initial values at B. Simi- 
larly, we can interpolate from B and the second node downstream of 
the shock (Q in Fig. 3) to get final values at the first node 
downstream (P in Fig. 3). 

The shock calculation is now completed. At the predictor 

XXX X 

level, however, we must re-evaluate fj , f 3 , f 5 and f £ at P, 
because these quantities are needed at the corrector level for 
y 

point Q (f,, can be left unchanged, since only g and h, as defined 
by (31), are used). The re-evaluation is possible using the actu- 
al increments of a, u and s in time and (23), (25), (26) and (24) 
XXX X 

to determine f 6 , fj , f s and f 3 , in that order. 


6.2. Formation and tracking of shocks . 

We opted for not attempting any tracking of oblique shocks with 
transitions from supersonic to supersonic flow. In the present class 
of problems, such shocks have a normal Mach number very close to 1 . 
Therefore, the entropy increase across the shock is minimal, and it 
can be neglected. Moreover, the shock is not influenced by down- 
stream conditions. Therefore, we may assume that it finds its loca- 
tion correctly. In conclusion, only the formation of shocks with 
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transition from supersonic to subsonic flow is to be monitored. 


Testing on shock formation begins only when there is at least 

X 

one grid cell across which changes from positive to negative. 
Each Y=constant line is swept from left to right, searching for a 

X 

possible cell where A . changes from positive to negative (and which 
has not been marked as a shock already). The node to the left of the 
cell is marked as a possible A-point of a shock. An approximate ini- 

X 

tial location of the shock is assumed in the cell, where \ 2 =0. New 
shock points are not generated if old shock points already exist 
along the same Y=constant line within two cells on either side. Once 
a shock point has been generated, the testing proceeds along the same 
Y=constant line until the entire line has been monitored. 

All shock points so generated over the computational grid are 
stored in a one- dimensional array. The next step consists of group- 
ing the shock points in successive links, so that each link contains 
shock points arranged in the order of increasing Y and the distance 
between two successive points of the same link is not greater than 
one cell. 

Once there is more than one point in one link, the direction of 
the normal at a point is easily found. At a wall, the normal to the 
shock is forced to be parallel to the wall. Otherwise, the normal is 
evaluated by using centered differences along the shock. 

If a shock point moves to an adjacent cell downstream, the 
values upstream are extended to the new upstream grid node; similar- 
ly, if the shock point moves to an adjacent cell upstream, the new 
downstream grid node is given the most recent values at the B-point. 
Such updatings are made at the end of the corrector level only; con- 
sequently, there is no need for re-evaluation of f-terms. 


6.3. Accuracy of the shock calculation . 

There is no doubt that the present technique of shock fitting is 
both fast and very accurate. Speed is achieved because the calcula- 
tion at any node remains the same as for a shockless flow, plus a 
minimal amount of additional work, as described in Section 6.1, per- 
formed at the shock points only. Even more important is the accuracy 
of the results. 

The entire calculation of the shock depends on two parameters, 
the shock angle and I. The angle is conveniently defined by centered 
differences, based on a shock geometry which is clearly described, 
regardless of the number of shocks and their relative positions. The 
other parameter, Z, as defined by (45), contain two variables comput- 
ed at point A, and one computed at point B. Let us consider f irst 
the case of Fig. 4, where three consecutive shock points remain 
within the same grid column. Since the flow is supersonic at A, all 

X Y 

the f -s are computed from upstream, correctly. All the f -s are 
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Figure 4.- Nearly vertical shock wave. 


computed as in any ordinary point, correctly. Therefore, the values 
of u^ and a A are correct. If we single out a g and Ug, which concur 

X 

to form R^g, we see that neither a nor u are correct, because they 
X 

depend on f 3 which is improperly computed across the shock; 

x 

nevertheless, the combination R 2B is correct, because the contribu- 

X XXX 

tion of fj cancels out (the contribution of f„ +f £ -2f 3 is minimal 

and can be set equal to zero). 



i i * i r 


Figure 5.- Shock waves passing through X=constant lines. 

Y 

In the case of Fig. 5, the terra, t 1 is evaluated at C using a 

Y Y Y 

local value of \ 1 and R 1 at Dj instead of at D, and the term, f z 

Y Y 

is evaluated at D using a local value of A. and R 2 at Cj instead of 
at C. Such devices can be loosely interpreted by saying that the 
derivatives are taken after subtracting the jumps across the shock. 
It is interesting to note that none of these devices requires the use 
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of IF-statements. 


7. A numerical example. 

From the several cases analyzed to test the present technique, 
we chose one which has some interesting features and the results of 
which can be compared to calculations by other authors. The geometry 
of the duct is shown in Fig. 1. The values of x 1 , x it x 0 and b c are 

x, = 1 .5 
x z = 2.5 
x c = 0.66937 
b 0 = 1.26775 

Therefore, the wedge has a 20$ slope. The Mach number at the entry 
section is 1.6. The calculation has been made starting from rest im- 
pulsively (this corresponds to the third model of Section 2, with an 
infinite acceleration). Two runs have been made, one with a 50 by 20 
mesh, and another with a 100 by 30 mesh. In both runs, the CFL 

number has been taken as 1 . In these runs, lengths have been scaled 
to the exit width, and speeds to the speed of sound at the entrance, 
divided by /T7 Therefore, times are scaled to the ratio between the 
width and the speed of sound at the entrance, multiplied by /Y?b 0 . 

The figures presenting results for these runs are self- 
explanatory. In Fig. 6, four sets of isomachs are shown, at increas- 
ing times, together with the shape of the shocks which have been fit- 
ted. Note the appearance of a second shock of very low strength, in 
the latter stages of the flow evolution. Note also the good resolu- 
tion and precision of details, despite the use of a mesh which is 
considered rather coarse in the existing literature. Better results, 
but not so much better as to justify the effort, are obtained by us- 
ing a finer mesh, as shown in the second set of plots, Fig. 7. In 
Fig. 8, a typical plot of entropy is shown, at the same step as the 
final plot of Fig. 7. Finally, in Fig. 9 we reproduce some figures 
from Ref. 4. Such figures have been obtained using a 150 by 50 mesh. 
We believe that the comparison shows that the present technique is 
more efficient. 

The code used so far has not been optimized; no attempt has been 
made yet to fully exploit the advantages of vectori zation. The com- 
putational time required on a CRAY 1 machine is of the order of 20 
seconds for 400 steps, using a 50 by 20 mesh. 
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Figure 6.- Isomachs at four times for 50 by 20 mesh. 
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Figure 7.- Isomachs at four times for 100 by 30 mesh. 
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